PROGRAM GAMMA1
!***********************************************************************************
!	This program uses the sigma profiles of two pure components to calculate 
!   the liquid-phase activity coefficients in a solution.  This is the first 
!   step in predicting VLE for mixtures.  
!
!	This program uses the COSMO-SAC model as published (Lin, S.T., 
!	S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913).
!
!	THIS PROGRAM WRITTEN BY:
!	  RICHARD OLDLAND (roldland@vt.edu)    MIKE ZWOLAK (zwolak@caltech.edu)
!	  DEPARTMENT OF CHEMICAL ENGINEERING   PHYSICS DEPARTMENT
!	  VIRGINIA TECH                        CALIFORNIA INSTITUTE OF TECHNOLOGY
!	  BLACKSBURG, VA 24060                 PASADENA, CA 91125
!
!	PHYSICAL CONSTANTS AND PARAMETERS:
!	  EO = PERMITTIVITY IN A VACUUM (e**2*mol/Kcal*Angstrom)
!	  AEFFPRIME = EFFECTIVE SURFACE AREA (ANGSTROMS**2) --FROM LIN
!	  RGAS = IDEAL GAS CONSTANT (Kcal/mol*K)
!	  VNORM = VOLUME NORMALIZATION CONSTANT (A**3) --FROM LIN
!	  ANORM = AREA NORMALIZATION CONSTANT (A**2) --FROM LIN
!	  COORD = THE COORIDINATION NUMBER --FROM LIN 
!	  CHB = HYDROGEN BONDING COEFFICIENT (Kcal/mole*Angstroms**4/e**2)
!	  SIGMAHB = CUTOFF VALUE FOR HYDROGEN BONDING (e/Angstrom**2)
!	  EPS = RELATIVE PERMITTIVITY --FROM LIN
!	  ALPHAPRIME = A CONSTANT USED IN THE MISFIT ENERGY CALCULATION
!
!	INPUT PARAMETERS:
!	  SYSTEMP = THE SYSTEM TEMPERATURE (K)
!	  COMP = NUMBER OF COMPONENTS IN THE SYSTEM --SET TO 2 FOR BINARY
!	  SYSCOMP = NAMES OF COMPONENTS IN THE SYSTEM
!	  VCOSMO = CAVITY VOLUME FROM COSMO OUTPUT (A**3)
!	  ACOSMO = MOLECULAR SURFACE AREA FROM COSMO OUTPUT (A**2) --THE SUM
!		OF THE INDIVIDUAL PROFILE.  
!
!
!	LITERATURE CITED:
!  Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the
!       Quantitative Calculation of Solvation Phenomena. J. Phys. Chem 1995, 99, 2224.
!  Klamt, A.; Jonas, V.; Burger, T.; Lohrenz, J. Refinement and Parameterization of 
!	COSMO-RS. J. Phys. Chem A 1998, 102, 5074.
!  Klamt, A.; Eckert, F.; COSMO-RS: A Novel and Efficient Method for the a Priori 
!	Prediction of Thermophysical Data of Liquids.  Fluid Phase Equilibria 2000, 
!	172, 43.
!  Lin, S.T.; Sandler, S. A Priori Phase Equilibrium Prediction from a Segment 
!       Contribution Solvation Model. Ind. Eng. Chem. Res, 2002, 41, 899 
!  Lin, S.T.;  Quantum Mechanical Approaches to the Prediction of Phase Equilibria: 
!	Solvation Thermodynamics and Group Contribution Methods, PhD. Dissertation, 
!	University of Delaware, Newark, DE, 2000
!
!
!	PROGRAM CURRENTLY SETUP FOR BINARY MIXTURES ONLY 
!**********************************************************************************
IMPLICIT NONE
REAL,PARAMETER:: EO = 2.395*10.0**-4, AEFFPRIME = 7.5, RGAS = 0.001987
REAL,PARAMETER:: VNORM = 66.69, ANORM = 79.53
REAL :: FPOL, ALPHA, ALPHAPRIME, COORD, EPS, SYSTEMP, SIGMAHB, CHB, FRAC1, FRAC2
REAL :: SYSPRES, SIGMAACC, SIGMADON, SUMMATION, BOTTHETA, BOTPHI, PHI1, PHI2
REAL :: THETA1, THETA2, L1, L2, GAMMASG1, GAMMASG2, GAMMA, GAMMA2, SUMGAMMA1
REAL :: SUMGAMMA2, N1, N2, LNGAMMA, LNGAMMA2
INTEGER :: I, J, K, L, M, COMPSEG, COUNT, COMP
REAL, DIMENSION(2):: VCOSMO, ACOSMO, RNORM, QNORM 
REAL, DIMENSION(:), ALLOCATABLE:: COUNTER, DENOM, PROFILE, NUMER, SEGGAMMA
REAL, DIMENSION(:), ALLOCATABLE:: SEGGAMMAOLD, CONVERG
REAL, DIMENSION(:,:), ALLOCATABLE:: SIGMA, DELTAW, SEGGAMMAPR, SEGGAMMAOLDPR, CONPR
CHARACTER (16), DIMENSION(2):: SYSCOMP
CHARACTER (256), DIMENSION(2):: FILENAME
CHARACTER :: ANSWER
CHARACTER (256):: OUTPUT

COMPSEG = 51 !NUMBER OF INTERVALS FOR THE SIGMA PROFILE
EPS = 3.667 !(LIN AND SANDLER USE A CONSTANT FPOL WHICH YEILDS EPS=3.68)
COORD = 10.0 !(KLAMT USED 7.2)  
SIGMAHB = 0.0084
CHB = 85580.0
COMP =2 

FPOL = (EPS-1.0)/(EPS+0.5)
ALPHA = (0.3*AEFFPRIME**(1.5))/(EO)
ALPHAPRIME = FPOL*ALPHA


ALLOCATE(SIGMA(COMPSEG,COMP), COUNTER(COMPSEG), PROFILE(COMPSEG), NUMER(COMPSEG),&
	DENOM(COMPSEG), DELTAW(COMPSEG,COMPSEG), SEGGAMMA(COMPSEG), SEGGAMMAOLD(COMPSEG), &
	CONVERG(COMPSEG), SEGGAMMAPR(COMPSEG,COMP), SEGGAMMAOLDPR(COMPSEG,COMP), &
	CONPR(COMPSEG,COMP))


!DEFINE SYSTEM TEMPERATURE (K)
!WRITE(*,*) "ENTER IN THE SYSTEM TEMPERATURE (K)"
!READ(*,*) SYSTEMP 
SYSTEMP=330.15

!DEFINE THE SYSTEM AS WELL AS THE AREA AND VOL FROM THE COSMO CALCULATION
!DO I=1, COMP
!	WRITE(*,*) "ENTER THE NAME OF COMPONENT", I
!	READ(*,*) SYSCOMP(I)
!	WRITE(*,*) "ENTER THE LOCATION AND NAME OF THE  ", SYSCOMP(I), "SIGMA PROFILE (256 character max)"
!	READ(*,*) FILENAME(I)
!	WRITE(*,*) "ENTER THE CAVITY VOLUME FOR COMPONENT",I
!	READ(*,*) VCOSMO(I)
!END DO
SYSCOMP(1)="methyl acetate"
FILENAME(1)="VT2005/VT2005-0638-PROF.txt"
VCOSMO(1)=97.00036
SYSCOMP(2)="water"
FILENAME(2)="VT2005/VT2005-1076-PROF.txt"
VCOSMO(2)=25.73454

!OPEN THE SIGMA PROFILES FOR THE PURE COMPONENTS
OPEN(UNIT=1, FILE=FILENAME(1), STATUS="OLD", ACTION="READ", POSITION="REWIND")
OPEN(UNIT=2, FILE=FILENAME(2), STATUS="OLD", ACTION="READ", POSITION="REWIND")

!READ INDIVIDUAL SIGMA PROFILES; COUNTER IS THE SIGMA VALUE, SIGMA IS P(SIGMA)
DO K = 1, COMP
	ACOSMO(K) =0.0
	DO J=1, COMPSEG
		READ(K,*)COUNTER(J), SIGMA(J,K)
		ACOSMO(K)=ACOSMO(K)+SIGMA(J,K)
	END DO
END DO
CLOSE(1)
CLOSE(2)

!ESTABLISH OUTPUT FILES.  THE PROGRAM CREATES A FILE FOR THE GAMMA-X DATA
WRITE(*,*) "ENTER THE LOCATION AND NAME FOR THE OUTPUT FILE, INCLUDE EXTENSION."
READ (*,*) OUTPUT
OPEN(UNIT=11, FILE = OUTPUT, STATUS="NEW")
WRITE(11,*) "SYSTEMP", SYSTEMP, "KELVIN"
5  FORMAT (1X,A10,5X,A12,5X,A12,5X,A12,5X,A12)
6  FORMAT (1X,A10,5X,A9,5X,A10,7X,A10,5X,A10)
WRITE(11,6) "MOLE FRAC", "GAMMA1", "GAMMA2", "LNGAMMA1", "LNGAMMA2"
WRITE(11,5) "X1", SYSCOMP(1), SYSCOMP(2), SYSCOMP(1), SYSCOMP(2)



!VARYING MOLE FRACTIONS // ONLY WORKS FOR BINARY MIXTURE
DO FRAC1 = 0.005, 0.995, 0.01
FRAC2 = 1.0 - FRAC1
!CALCULATE THE MIXTURE SIGMA PROFILE 
	DO J =1,COMPSEG
		NUMER(J) = FRAC1*SIGMA(J,1) + FRAC2*SIGMA(J,2)
		DENOM(J) = FRAC1*ACOSMO(1) + FRAC2*ACOSMO(2)
		PROFILE(J) = NUMER(J)/DENOM(J)

	END DO
	DO I = 1, COMPSEG
		DO K = 1, COMPSEG
			IF (COUNTER(I)>=COUNTER(K)) THEN
				SIGMAACC = COUNTER(I)
				SIGMADON = COUNTER(K)
			END IF
			IF (COUNTER(I)<COUNTER(K)) THEN
				SIGMADON = COUNTER(I)
				SIGMAACC = COUNTER(K)
			END IF
			
		DELTAW(I,K) = (ALPHAPRIME/2.0)*(COUNTER(I)+COUNTER(K))**2.0 + CHB *   &
			MAX(0.0,(SIGMAACC - SIGMAHB))*MIN(0.0,(SIGMADON + SIGMAHB))
		END DO
	END DO


!ITERATION FOR SEGMENT ACTIVITY COEF. (MIXTURE)
SEGGAMMA = 1.0
L = 0
DO 
	SEGGAMMAOLD = SEGGAMMA
	DO I = 1, COMPSEG
		SUMMATION = 0.0
		DO K = 1, COMPSEG
			SUMMATION = SUMMATION + PROFILE(K)* SEGGAMMAOLD(K) * &
				EXP(-DELTAW(I,K)/(RGAS*SYSTEMP))
		END DO
		SEGGAMMA(I)=EXP(-LOG(SUMMATION))
!		SEGGAMMA(I)=1.0/SUMMATION
		SEGGAMMA(I)=(SEGGAMMA(I)+SEGGAMMAOLD(I))/2.0
	END DO
	DO I=1, COMPSEG
		CONVERG(I)=ABS((SEGGAMMA(I)-SEGGAMMAOLD(I))/SEGGAMMAOLD(I))
	END DO
L = L+1
IF (MAXVAL(CONVERG) <=0.000001) EXIT
END DO
WRITE(*,*) "SEGGAMMA niter:", L


!ITERATION FOR SEGMENT ACITIVITY COEF (PURE SPECIES)
DO L = 1, COMP
	SEGGAMMAPR (:,L) = 1.0
	DO 
		SEGGAMMAOLDPR (:,L) = SEGGAMMAPR (:,L)
		DO I = 1, COMPSEG
			SUMMATION = 0.0
			DO K = 1, COMPSEG
				SUMMATION = SUMMATION + (SIGMA(K,L)/ACOSMO(L))*SEGGAMMAOLDPR(K,L) * &
					EXP(-DELTAW(I,K)/(RGAS*SYSTEMP))
			END DO
			SEGGAMMAPR(I,L)=EXP(-LOG(SUMMATION))
			SEGGAMMAPR(I,L)=(SEGGAMMAPR(I,L)+SEGGAMMAOLDPR(I,L))/2.0
		END DO
		DO I=1, COMPSEG
			CONPR(I,L)=ABS((SEGGAMMAPR(I,L)-SEGGAMMAOLDPR(I,L))/SEGGAMMAOLDPR(I,L))
		END DO
		
	IF (MAXVAL(CONPR) <=0.000001) EXIT
	END DO
END DO

!THE STAVERMAN-GUGGENHEIM EQUATION
DO I = 1,COMP
	RNORM(I) = VCOSMO(I)/VNORM
	QNORM(I) = ACOSMO(I)/ANORM
END DO

BOTTHETA = FRAC1*QNORM(1) + FRAC2*QNORM(2)
BOTPHI = FRAC1*RNORM(1) + FRAC2*RNORM(2)

THETA1 = (FRAC1*QNORM(1))/BOTTHETA
THETA2 = (FRAC2*QNORM(2))/BOTTHETA

PHI1 = (FRAC1*RNORM(1))/BOTPHI
PHI2 = (FRAC2*RNORM(2))/BOTPHI

L1 = (COORD/2.0)*(RNORM(1)-QNORM(1))-(RNORM(1)-1.0)
L2 = (COORD/2.0)*(RNORM(2)-QNORM(2))-(RNORM(2)-1.0)

!GAMMASG1 AND GAMMASG2 ARE ACTUALLY LNGAMMASG
GAMMASG1 = LOG(PHI1/FRAC1)+(COORD/2)*QNORM(1)*LOG(THETA1/PHI1)+L1-(PHI1/FRAC1)* &
	(FRAC1*L1 + FRAC2*L2)
GAMMASG2 = LOG(PHI2/FRAC2)+(COORD/2)*QNORM(2)*LOG(THETA2/PHI2)+L2-(PHI2/FRAC2)* &
	(FRAC1*L1 + FRAC2*L2)


!CALCULATION OF GAMMAS
SUMGAMMA1 = 0.0
SUMGAMMA2 = 0.0
DO I = 1, COMPSEG
	SUMGAMMA1 = SUMGAMMA1 +((SIGMA(I,1)/AEFFPRIME)*(LOG(SEGGAMMA(I)/(SEGGAMMAPR(I,1)))))
	SUMGAMMA2 = SUMGAMMA2 +((SIGMA(I,2)/AEFFPRIME)*(LOG(SEGGAMMA(I)/(SEGGAMMAPR(I,2)))))
END DO


GAMMA =EXP(SUMGAMMA1 + (GAMMASG1))
GAMMA2=EXP(SUMGAMMA2 + (GAMMASG2))
LNGAMMA = LOG(GAMMA)
LNGAMMA2 = LOG(GAMMA2)


WRITE(11,*) FRAC1,GAMMA, GAMMA2, LNGAMMA, LNGAMMA2


END DO


END PROGRAM GAMMA1
  
